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1.  INTRODUCTION 


Increasingly  there  is  a  need  to  more  accurately  predict  the  viscous 
resistance  on  ship  hulls,  particularly  in  light  of  modern-day  fuel 
costs.  The  rapid  increase  in  fuel  costs  cries  out  as  a  mandate 
to  ship-hull  designers  to  reduce  hull  resistance,  both  wave-making 
and  viscous.  Experience  has  shown  that  significant  advances  in 
design  efficency  often  require  more  sophisticated  design  tools, 
and  in  today's  environment,  particularly  those  of  the  computational 
variety . 

Airplane  and  missile  designers  have  had  a  great  deal  of  success  in 
developing  and  utilizing  the  concept  of  the  "numerical  wind  tunnel". 
As  both  computational  algorithms  and  accurate  engineering  sets  of 
equations  describing  turbulent  fluctuations  and  stresses  have  been 
developed  and  improved,  so  has  the  ability  to  replace  expensive 
wind-tunnel  tests  by  a  series  of  computer  runs.  In  addition  to 
being  relatively  inexpensive,  such  numerical  simulations  can  more 
exactly  simulate  full-scale  flow  conditions  than  can  wind-tunnel 
tests . 

In  the  ship  context  there  is  just  as  pressing  a  need  to  replace 
expensive  tests  by  numerical  simulations.  In  fact,  in  light  of  the 
inherent  contradiction  involved  in  trying  to  simultaneously  match 
Froude  and  Reynolds  numbers,  then  is  an  even  more  pressing  need  to 
utilize  the  "numerical  towing  tank"  in  the  design  process. 

The  state  of  the  art  of  numerical  simulations,  including  viscous 
effects,  is  further  advanced  in  the  aerodynamic  context  then  it  is 
in  the  hydrodynamic  case.  Conceptually  this  is  sensible  as  the 
flow  past  a  ship  hull  involves  a  free  surface  and  is  inherently 
three  dimensional  while,  by  contrast,  much  of  the  flow  area  over 
an  airplane  wing  involves  a  homogeneous  fluid  and  can  be  treated 
as  two  dimensional.  This  two  dimensionality  has  been  of  especial 
value  to  aerodynamicists  as  they  have  been  able  to  justify  focusing 
on  two-dimensional  flows.  A  great  deal  of  progress  in  Improving 
computational  efficiency  has  been  made  because  of  the  relative 


simplicity  of  a  two-dimensional  computation,  progress  which  would 

have  been  much  harder  to  achieve  if  the  aerodynamicist ’ s  primary 

application  had  demanded  focusing  on  three  dimensionality.  As  a 

result,  computing  separated  flows  can  be  done  today1  at  less  than 

2 

a  percent  of  the  cost  seven  or  eight  years  ago  .  And,  of  course, 
these  methods  are  now  being  implemented  in  three  dimensions. 

Even  in  three  dimensions  the  aerodynamicist  has  an  easier  time  of 
it  than  the  hydrodynamicist  because  he  needn’t  deal  with  a  free 
surface.  Nevertheless,  the  day  of  the  numerical  towing  tank  draws 
nearer  and  nearer.  The  purpose  of  this  project  has  been  to  help 
in  hastening  the  arrival  of  the  numerical  towing  tank. 

To  accomplish  this  end,  we  have  addressed  the  viscous  portion  of 
the  problem.  We  have  applied  a  widely-tested  two-equation  model 
of  turbulence  (Section  2)  to  (a)  real  ship  hulls  and  (b)  "thick" 
turbulent  boundary  layers  on  submerged  axisymmetric  bodies.  The 
object  of  the  ship-hull  computations  has  been  to  test  the  model's 
accuracy  in  the  limit  of  classical  "thin-shear-layer"  approximations. 
The  latter  computations  have  been  done  as  an  attempt  to  pinpoint 
the  special  characteristics  of  and  to  test  the  model  for  "thick" 
boundary  layers. 

In  performing  the  analysis  we  have  developed  a  new  three-dimensional 
boundary-layer  program  appropriate  to  arbitrary  ship  hulls.  A  great 
deal  of  the  effort  in  this  project  focused  on  devising  an  accurate 
numerical  procedure  compatible  with  the  turbulence  model  equations. 
Section  3  presents  details  of  the  algorithm  including  a  detailed 
accuracy  study. 

In  Section  4  we  present  results  of  the  two  ship-hull  computations, 
including  comparison  of  computed  and  measured  flow  properties. 

Section  5  includes  details  of  our  thick-boundary-layer  study.  Re¬ 
sults  and  conclusions  follow  in  Section  6.  The  Appendix  includes 
computer-program  listings  and  an  explanation  of  program  input  and 
output . 


2.  EQUATIONS  OF  MOTION 

■3 

In  this  section  we  describe  the  Wilcox-RubesinJ  two-equation  model 
of  turbulence.  First  we  summarize  the  equations  of  motion  in¬ 
cluding  all  closure  coefficients.  Then,  we  state  the  boundary 
conditions  appropriate  to  smooth  and  rough  surfaces. 

2.1  THE  TURBULENCE  MODEL 

For  an  incompressible  fluid  of  density,  p,  and  kinematic  viscosity, 
v,  we  must  solve  the  equations  of  mass  and  momentum  conservation, 
viz, 

3u. 

TT~  *  0  (1 


.  _  i  3£  _s_  r 

uj  axj  p  3x^  {_  v  3xj 


where  steady  flow  conditions  have  been  assumed.  In  Equations  (1-2), 
ui  and  denote  velocity  and  position  vector,  p  is  static  pressure 
and  is  the  Reynolds-stress  tensor.  To  close  the  system,  is 

assumed  proportional  to  the  mean  strain  rate  tensor,  ,  wherefore 

Tij  *  2eSiJ  '  fk6ij  (3) 

where  3u,  3u. 


sij  s  15  +  33^}  1 

The  quantity  k  is  turbulent  mixing  energy  and  e  is  kinematic  eddy 
viscosity.  The  latter  is  defined  in  terms  of  k  and  the  turbulent 
dissipation  rate,  w,  by 

e  =  Y*k/w  i 

where  y*  is  a  closure  coefficient.  Finally,  the  turbulence  para¬ 
meters  k  and  w  satisfy  the  following  equations. 


uj%-  hj  nj-  8*“k  +  %[<u  +  °*e)  %] 


-3- 


.  2  2  3U 

du  U)  _ 

uj3Xj  Y  k  Tij  3x 


J-[6  +  20  <%>’]  +  Hj[<''+0£)  lij] 


where  i  is  the  turbulent  length  scale  defined  by: 

l  =  ( 8 ) 

Equations  (5-7)  contain  six  closure  coefficients,  viz,  B,  3*,  Y,  Y * , 
o  and  c*.  The  values  of  these  coefficients  have  been  determined 
from  very  general  arguments  based  on  widely  observed  properties  of 
turbulent  flows.  The  values  used  in  the  present  study  are: 


B  =  3/20  ;  B*  = 

a  =  h  >  o*  = 

Y*  =  1  -  ( 1—  X 2 )exp(-ReT) 

YY*  =  [l-(l-X2  )exp(-ReT/^)] 

where  ^  =  lT 


9/100 

h 


In  Equations  (9),  the  quantity  Re,p  is  the  Reynolds  number  of  the 
turbulence  defined  by 


2 . 2  BOUNDARY  CONDITIONS 


Application  of  the  no-slip  boundary  condition  yields 


u^  =  0 


k  =  0 


at  solid  boundaries 


(11) 


The  boundary  conditions  for  w  appropriate  to  smooth  and  rough  sur- 
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faces  have  been  shown  by  Wilcox  and  Traci  to  be 


as  y  -*■  0 


Smooth  Surface 


S  ut 


at  y  =  0  ;  Rough  Surface 


where  uT  is  friction  velocity  and  S  is  a  function  of  dimensionless 


roughness  height  k  =  uTk/v  defined  by 

S  =  (36/k+)2  +  (8/k+)?5 

Equations  (1-13)  completely  define  the  turbulence  model  and 
surface  boundary  conditions.  All  analysis  in  this  project  ha 
been  done  using  Equations  (1-13). 


3.  NUMERICAL  CONSIDERATIONS 


Because  our  primary  interest  in  this  project  is  in  three-dimensional 
boundary  layers,  a  substantial  effort  has  gone  into  devising  an 
accurate,  stable  and  efficient  numerical  procedure  for  solving 
Equations  (1-13)  in  the  3-D  boundary  layer.  Few  computer  programs 
have  been  developed  with  sufficient  generality  for  arbitrary  ship- 
hull  applications  and  we  have,  in  effect,  developed  such  a  pro¬ 
gram  from  scratch. 

Cur  original  starting  point  was  a  generalized  version  of  the  3-D 

5 

boundary-layer  program  developed  by  Cebeci,  et  al  .  The  r.  >rical 
procedure  embodied  in  the  program  is  based  on  Keller' s^  "  x  Method" 
After  a  great  deal  of  numerical  experimentation  and  frust  tion, 
we  have  concluded  that  the  Box  Method  is  incompatible  wit  >  ? 

turbulence-model  equations. 

This  is  not  to  say  the  Box  Method  is  deficient  as  a  numerical  scheme 
On  the  contrary,  the  method  has  proven  to  be  very  efficient  and  very 
accurate  in  many  applications.  However,  advanced  turbulence  models 
such  as  the  one  being  used  often  require  special  numerical  treatment 
Over  the  years,  we  have  found  that  many  excellent,  proven,  numerical 
schemes  just  don't  work  well  in  solving  the  equations  attending 
such  a  model.  We  now  know  the  Box  Method  falls  into  this  category. 

To  remedy  this  problem,  we  have  completely  rewritten  the  program. 

It  now  embodies  a  numerical  method  which  has  proven  compatible  with 
the  turbulence  model.  In  the  following  sections  we  first  discuss 
the  new  numerical  method.  Next,  we  show  typical  numerical  results 
for  a  flat-plate  boundary  layer,  including  comparisons  with  exper¬ 
imental  data.  Finally,  we  present  results  of  a  numerical  accuracy 
study. 

3.1  THE  NUMERICAL  PROCEDURE 

For  many  years  we  have  obtained  accurate  two-dimensional  boundary 
layer  solutions  for  our  model  equations  using  the  classical  Blot- 
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tner  method.  Hence,  we  decided  to  extend  the  method  to  three 


-6- 


dimensions  in  hopes  it  wo"1  i 


retain  its  desirable  stability  and 
accuracy  cho^„-terlstics .  Before  embarking  on  such  a  major  re¬ 
programming  effort  we  made  a  brief  review  of  Blottner's  more  re¬ 
cent  work.  We  found  that,  as  an  improvement  to  his  original  pro- 

g 

cedure ,  Blottner0  has  devised  his  "Variable  Grid"  algorithm.  This 
revised  method  offers  far  greater  accuracy  than  his  original  method, 
particularly  for  very  coarse  finite-difference  grids.  We  decided 
to  use  the  improved  algorithm. 


In  essence,  the  "Variable  Grid"  method  uses  a  conservation  form 
treatment  for  diffusion  terms.  That  is,  in  the  original  Blottner 
method,  diffusion  terms  must  first  be  expanded  according  to: 


9  /  9u  N  _  9  2u  9_e_  du 
9y  3y '  %y2  9y  9  y 


(1^4 


Then,  central  difference  approximations  are  used  for  9u/9y,  9e/9y, 
and  92u/9y2.  By  contrast,  in  the  improved  Blottner  method,  we  write 
directly  (subscripts  denote  mesh  point  number): 

ri_(E3u)]  .  Su3 *h  -  h-’i  ‘“j-'i  (15 

L9y^  9yjJj  V(Ay .  +  Ay  ._1) 


6u 


J+*S 


(u 


J+l  -  Uj)/Ayj 


(1 


and  similarly  for  (<5u)j_^.  The  quantity  Ay^  =  yj  +  1  -  yj  is  the 
incremental  change  in  the  normal  coordinate  y.  As  will  be  shown 
in  Subsection  3-3,  this  straightforward  modification  permits  much 
greater  stretching  of  the  grid  normal  to  the  surface  than  is  pos¬ 
sible  with  the  original  Blottner  method. 


Thus  far,  we  have  spoken  only  of  the  direction  normal  to  the  sur¬ 
face.  Because  the  three-dimensional  boundary-layer  equations  have 
a  real  characteristic  which  is  not  necessarily  aligned  with  the 
freestream  flow  direction,  any  accurate  numerical  procedure  must 
be  consistent  with  the  attending  zone  of  dependence.  The  original 
Box  Method  does  not  treat  this  problem  properly  and  to  avoid 


future  difficulties  we  decided  to  accomodate  the  zone-of-dependence 
feature  in  the  revised  numerical  scheme.  In  an  unpublished  study, 

Q 

Wilcox  and  James  devised  a  family  of  explicit,  unconditionally 
stable  marching  methods  which  potentially  could  deal  with  this 
problem.  A  review  of  that  study  showed  that  most  of  the  family 
were  of  first-order  accuracy.  However,  one  of  the  most  promising 
of  the  family  is  second-order  accurate.  Another  brief  review  of 
boundary-layer  literature  showed  that  this  scheme  has  been  used  by 
many  others  in  three-dimensional  boundary-layer  applications.  The 
scheme  is  attributed  to  Krause1  .  We  have  opted  to  use  Krause 
scheme  in  the  revised  program. 


The  essence  of  the  Krause  procedure  is  the  manner  in  which  cross- 
flow  convective  terms  are  treated.  As  an  illustration,  the  term 
w3u/3z  (where  u  and  w  denote  streamwise  and  crossflow  velocity 
components  and  z  is  the  crossflow  coordinate)  is  discretized 
according  to: 


u, 


,  „  \  n  /“i-l.k+l  "  ui-l ,k  .  ui,k  “  ui,k-l 

(w9u/3j>1Jljk  -  <s  - *-  +  -  - 


Jk-1 


(1 


where  i  and  k  denote  mesh  point  number  in  the  x  and  z  directions. 


respectively,  Azk  =  z 
is  as  shown  below. 


k+1 


-  z,  ,  and  the  finite-difference  molecule 
k  * 


(i-l,k+l) 


(i-1 ,k) 
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Centered  Here 
(i-H»k) 

(i  »k) 

(i.k-1) 


This  scheme  is  unconditionally  stable  for  positive  crossflow 
(w  >  0)  and  conditionally  stable  for  negative  crossflow  (w  <  0), 
the  stability  condition  being: 


Thus,  for  equally  spaced  mesh  points  (Ax  =  Az),  this  scheme  is 
stable  for  negative  crossflow  angles  up  to  ^5  degrees  from  the 
streamwise  direction. 

3.2  PRELIMINARY  TEST  OF  THE  REVISED  PROGRAM 

A  key  difficulty  encountered  with  the  Box  Method  was  the  presence 
of  large  round-off  errors  in  the  equation  for  w2  .  Even  when  we 
devised  a  special  procedure  to  eliminate  these  errors  (caused  by 
the  numerical  method's  inability  to  accurately  compute  the  rapid 
variation  in  w2  as  y  -*•  0),  nontrivial  oscillations  in  properties 
such  as  skin  friction  remained.  All  of  these  difficulties  were 
encountered  for  the  simplest  boundary-layer  of  all,  the  two- 
dimensional  flat-plate  (constant  pressure)  boundary  layer.  Hence, 
as  a  preliminary  test  of  our  revised  3-D  boundary-layer  program, 
which  we  will  refer  to  as  EDDY3,  we  focus  on  the  flat-plate  boundary 
layer.  We  have  done  two  somputations ,  the  first  for  a  turbulent 
case  and  the  second  for  a  transitional  case. 

The  first  computation  was  initiated  from  approximate  profiles  at 
a  plate-length  Reynolds  number,  Rex,  of  one  million.  Using  *40 
equally  spaced  streamwise  steps,  the  computation  proceeds  to  a 
plate-length  Reynolds  number  of  five  million.  This  means  the  ratio 
of  streamwise  stepsize  to  boundary-layer  thickness  ranges  from  1.3 
to  6.8.  Such  large  steps  are  comparable  to  those  used  in  mixing- 
length  computations.  Only  61  mesh  points  were  used  normal  to  the 
surface  with  the  grid-point  incremental  spacing  increasing  in  a 
geometric  progression  with  a  progression  ratio  of  r  =  1.11. 

Figure  1  compares  computed  skin  friction,  cf,  and  momentum  thick¬ 
ness  Reynolds  number,  Re0,  with  the  Karman-Schoenherr  correlation. 

As  shown,  both  c^.  and  Re0  virtually  duplicate  the  correlation. 
Figures  2  through  4  compare  computed  velocity,  turbulent  dissipa¬ 
tion  rate  and  turbulent  mixing  energy  profiles  with  corresponding 
exact  theoretical  asymptotic  sublayer  and  wall-layer  (log  region) 
profiles.  In  all  cases,  the  numerical  predictions  clearly  fall 
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Figure  2.  Comparison  of  computed  velocity  profile  for  a  flat-plate  boundary  layer  with 
theoretical  asymptotic  sublayer  and  wall-layer  profiles;  61  mesh  points  and 


Figure  3 


Comparison  of  computed  specific  dissipation  rate 
profile  for  a  flat-plate  boundary  layer  with 
theoretical  asymptotic  sublayer  and  wall-layer 
profiles;  6l  mesh  points  and  r  =  1.11. 
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Figure  4.  Comparison  of  computed  turbulent  mixing  energy  profile 

for  a  flat-plate  boundary  layer  with  theoretical  asymptotic 
sublayer  and  wall-layer  profiles;  6l  mesh  points  and  r=l.ll 
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very  close  to  the  asymptotic  profiles,  thus  verifying  the  accuracy 
of  the  revised  numerical  procedure  in  the  absence  of  crossflow. 

For  our  second  test  we  begin  with  laminar  flow  and  integrate  all 
the  way  through  transition  to  test  whether  cf  oscillations  attend 
such  a  computation.  As  shown  in  Figure  5,  the  revised  program 
displays  no  such  oscillations.  Again,  the  geometric  progression 
ratio  used  was  r  =  1.11. 

In  conclusion,  the  success  achieved  in  these  two  flat-plate  cases 
leaves  us  confident  that  we  have  a  numerical  procedure  which  is 
at  once  stable  and  accurate.  While  these  computations  fail  to 
test  the  program’s  ability  to  handle  reverse  crossflow,  the  ship- 
hull  applications  of  Section  4  indicate  no  difficulties  are  en¬ 
countered,  even  with  very  large  crossflow  velocities. 

2.3  NUMERICAL  ACCURACY  STUDY 

Before  proceeding  to  more  complex  applications,  it  is  instructive 
to  examine  the  sensitivity  of  EDDY3  computations  to  grid-point 
spacing.  One  of  the  most  remarkable  features  of  the  revised 
program  is  that  in  past  two-dimensional  computations,  at  least  200 
mesh  points  would  have  been  needed  to  achieve  mesh-independent  solu¬ 
tions  for  an  Rex  of  five  million.  This  is  true  because  the  original 
Blottner  method  is  of  questionable  accuracy  for  values  of  r  in  ex¬ 
cess  of  1.0-4.  To  test  the  program’s  accuracy,  we  performed  a  series 
of  turbulent  flat-plate  boundary-layer  computations  in  which  r  was 
increased  all  the  way  to  1.^0.  To  establish  a  baseline  solution, 
we  used  our  two-dimensional  program  to  compute  the  flow  of  Sub¬ 
section  3-2  with  200  mesh  points  and  with  r  =  1. 03- 

Figure  6  shows  the  percent  error  in  skin  friction  as  a  function  of 
the  progression  ratio  r .  Remarkably,  even  with  an  r  as  large  as 
1.40  the  cf  error  is  only  8.5  percent.  Values  of  r  less  than  1.14 
yield  skin  friction  errors  less  than  2  percent. 
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Figure  7  recasts  the  percent  error  in  skin  friction  in  terms  of 
mesh  point  number.  As  shown,  solutions  are  more  or  less  mesh 
independent  beyond  about  60  mesh  points. 


4.  SHIP-HULL  APPLICATIONS 


We  turn  now  to  computation  of  the  three-dimensional  boundary  layer- 
on  two  ship  hulls,  viz,  the  SSPA  Model  720  and  the  HSVA  Tanker. 
These  are  the  two  test  cases  included  in  the  1980  SSPA-ITTC 
Workshop  on  ship  boundary  layers^.  For  the  SSPA  Model  720  we 
perform  three  separate  computations  to  test  EDDY3's  ability  to 
handle  reverse  crossflow  and  sensitivity  for  grid  point  spacing. 
Then,  the  HSVA  Tanker  computation  is  done  just  once  on  an  optimized 
grid. 

The  surface  finite  difference  grid  is  nonorthogonal  and  PROGRAM 
SHPMSH  (see  Section  A.l  of  the  Appendix)  is  used  to  generate  the 
grid.  Also,  PROGRAM  VELOC  (see  Section  A. 2  of  the  Appendix)  is 
used  to  interpolate  inviscid  velocities  onto  the  nonorthogonal 
grid. 


4.1  SSPA  MODEL  720 

4.1.1  Effect  of  Girthwise  Integration  Direction 

All  of  the  computations  performed  use  the  Workshop-supplied  Douglas- 
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Neumann  inviscid  velocity  distribution.  In  all  cases,  computation 
is  Initiated  from  fully-turbulent  boundary-layer  profiles  at  2x/L  = 
-0.6  (L  Is  the  hull  half  length).  The  profiles  used  match  the 
measured  momentum  thickness  and  skin  friction.  For  the  computations 
of  this  Subsection,  the  finite-difference  mesh  consists  of  31 
equally-spaced  points  extending  from  2x/L  =  -0.6  to  2x/L  =  0.9, 

11  equally-spaced  points  (in  terms  of  girth)  in  the  girthwise  direc¬ 
tion,  and  an  average  of  40  points  normal  to  the  surface.  Thus,  the 
mesh  consists  of  about  22,000  points.  Consistent  with  the  Douglas- 
Neumann  computation,  both  the  keel  line  and  the  waterline  have 
been  treated  as  symmetry  planes. 

To  further  test  the  program's  numerical  formulation,  we  perform 
the  computation  in  two  different  ways.  First,  we  perform  the 
girthwise  integration  starting  at  the  keel  and  integrating  toward 
the  waterline  on  each  streamwise  plane  (section).  Then,  we  repeat 
the  computation  Integrating  from  waterline  to  keel.  Doing  both 
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computations  checks  numerical-algorithm  consistency  and  the  pro¬ 
gram's  ability  to  handle  reverse  crossflow. 

Inspection  of  key  flow  properties  shows  virtually  no  difference 
between  the  two  solutions  through  the  entire  flowfield.  This 
confirms  that  the  program  is  indeed  able  to  handle  reverse  cross- 
flow  in  a  stable  and  accurate  manner.  It  also  confirms  that  over¬ 
all  the  algorithm  is  consistent  witi.  the  parabolic  nature  of  the 
boundary-layer  equations. 

*4.1.2  Effect  of  Mesh  Refinement 

In  the  next  computation  we  use  81  mesh  points  in  the  streamwise 
direction,  half  of  which  lie  between  2x/L  =  0.5  and  2x/L  =  0.9. 

A  total  of  21  equally-spaced  points  lie  in  the  girthwise  direc¬ 
tion  with  an  average  of  75  points  normal  to  the  surface.  This 
mesh  consists  of  about  128,000  points,  almost  six  times  the 
number  used  in  the  computations  of  Subsection  4.1.1.  The  girth- 
wise  integration  is  carried  out  from  keel  to  waterline. 

A  detailed  comparison  of  corresponding  flow  properties  for  the 
22,000-  and  the  128,000-mesh-point  computations  shows  that  over¬ 
all  flow  properties  differ  by  about  3  percent.  Even  in  regions 
of  rapid  change,  differences  never  exceed  5  percent.  These  ob¬ 
servations  are  consistent  with  our  accuracy  study  for  the  flat- 
plate  boundary  layer  which  indicate  that,  on  the  one  hand,  using 
40  points  normal  to  the  surface  yields  a  little  over  3  percent 
error  while,  on  the  other  hand,  using  75  points  results  in  about 
0.6  percent  error.  The  5  percent  discrepancies  almost  certainly 
attend  the  difference  in  streamwise  resolution  between  the  two 
computations . 

Results  of  this  numerical  test  indicate  there  is  little  point  in 
using  a  mesh  as  fine  as  the  one  with  128,000  points  if  economy 
is  a  factor.  That  is,  our  22,000-point  computations  require  ap¬ 
proximately  20  minutes  on  a  UNIVAC  1108  while  the  128,000-point 
computation  requires  about  90  minutes.  Note  that  the  Increase 
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in  machine  time  is  slower  than  linear  because  use  of  the  finer 
mesh  is  attended  by  fewer  iterations  in  the  solution  procedure 
at  each  mesh  point. 

4.1.3  Comparison  of  Computed  and  Measured  Momentum  Thickness 

In  all  three  computations,  we  predict  boundary-layer  separation 
over  a  portion  of  the  hull  beyond  2x/L  =  0.77.  Although  no  sep¬ 
aration  appears  to  have  been  observed  experimentally,  one  should 
also  note  that  the  experimentally  observed  boundary  layers  did 
not  experience  as  strong  an  adverse  pressure  gradient  as  those 
predicted  by  the  Douglas-Neumann  computation. 

Figure  8  compares  computed  and  measured  momentum  thickness,  0, 
on  three  lines  along  the  hull.  Line  A  is  the  keel  line.  Line  E 
is  a  line  well  below  the  waterline  where  the  boundary  layer  is 
more-or-less  of  the  classical  "thin"  variety.  Line  C  is  closer 
to  the  waterline  and  the  boundary  layer  approaches  the  more  com¬ 
plicated  "thick"  structure. 

Along  Line  A,  the  boundary  layer  is  truly  three  dimensional  as 
exhibited  by  its  curious  behavior  approaching  the  stern.  Speci¬ 
fically,  despite  entering  a  region  of  adverse  pressure  gradient, 
the  momentum  thickness  decreases.  This  behavior  occurs  because 
of  large  flow  divergence  near  the  stern.  While  the  computed  0 
lies  about  25  percent  above  measured  values,  note  that  our  pre¬ 
dicted  curve  is  much  closer  to  the  data  than  those  of  virtually 
all  Computors  in  the  SSPA-ITTC  Workshop  on  Ship  Boundary  Layers, 
most  of  whose  curves  were  high  by  a  factor  of  two! 

Along  Line  B,  we  exhibit  somewhat  larger  differences  from  the 
measured  values.  Computors  in  the  Conference  generally  came 
closest  to  experimental  data  on  this  line.  Our  prediction  is 
about  average  relative  to  Conference  participants'. 
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On  Line  C,  our  predicted  curve  follows  measured  values  quite 
closely  up  to  about  2x/L  =  0.5;  beyond  this  point  we  predict  a 
slower-than-measured  increase  in  momentum  thickness.  This  is 
unsurprising  as  no  provision  has  been  made  to  accomodate  any 
"thick"  boundary-layer  phenomena.  Our  prediction  is  much  closer 
to  the  data  than  most  Conference  participants'. 

1* .  2  HSVA  TANKER 

As  with  the  SSPA  Model  720,  we  use  the  Workshop-supplied  Douglas- 
Neumann  inviscid  velocity  distribution.  Computation  is  initiated 
from  fully-turbulent  boundary-layer  profiles  at  2x/L  =  -0.80.  Again, 
the  profiles  used  match  the  measured  momentum  thickness  and  skin 
friction.  The  finite-dif ference  mesh  consists  of  70  points  in  the 
streamwise  direction  extending  from  2x/L  =  -0.80  to  2x/L  =  0.90, 

15  equally-spaced  points  (in  terms  of  girth)  in  the  girthwise 
direction,  and  an  average  of  60  points  normal  to  the  surface.  Figure 
9  compares  computed  and  measured  momentum  thickness  on  three  sections 
viz,  for  2x/L  »  -0.7^4,  0.291  and  0.5 02.  The  numerical  boundary 
layer  separated  over  a  substantial  portion  of  the  hull  at  2x/L  =  0.75 
so  that  comparison  with  momentum-thickness  data  at  2x/L  =  0.88^  is 
not  possible.  Note  that  the  abscissa,  z,  is  the  nondimensional 
girth  with  z=0  on  the  keel  and  z=l  on  the  waterline.  Overall 
differences  between  computed  and  measured  0  are  similar  to  those 
obtained  for  the  SSPA  Model  720. 

In  summary,  we  have  tested  EDDY3  on  two  ship  hulls  and  encountered 
no  significant  numerical  difficulty.  Separation  on  both  hulls  is 
almost  certainly  caused  by  the  use  of  the  Neumann  inviscid  velocities 
so  separation  is  discounted  as  a  numerical  problem.  The  differences 
observed  between  computed  and  measured  momentum  thickness,  although 
generally  much  smaller  than  those  of  Workshop  participants,  are 
too  large  for  general  engineering  applications.  Thus,  our  predictive 
ability  requires  further  improvement  before  a  ship  designer  can 
apply  a  program  such  as  EDDY3  with  any  degree  of  confidence. 


Inspection  of  Figures  8  and  9  indicates  our  predictions  differ 
from  measurements  most  significantly  as  we  approach  the  stern. 

This  is  unsurprising  because,  on  the  one  hand,  we  expect  the 
boundary  layer  to  become  "thick"  so  that  classical  thin-shear-layer 
approximations  become  suspect  while,  on  the  other  hand,  our  compu¬ 
tations  have  been  done  with  classical  thin-shear-layer  approxima¬ 
tions.  In  the  next  section,  we  address  the  question  of  how  well 
our  model  equations  predict  properties  of  "thick"  boundary  layers 
if  we  abandon  the  thin-shear-layer  limit. 


5.  ANALYSIS  OF  THE  THICK  BOUNDARY  LAYER 


In  this  section  we  address  the  thick  boundary  layer  by  computing 
three  axisymmetric  flows  in  which  boundary-layer  thickness  becomes 
large  relative  to  the  body  radius.  First  we  use  perturbation  methods 
to  show  that  model  equation  solutions  are  consistant  with  the  Rao 
and  Richmond  scaling  laws  for  axisymmetric  boundary  layers.  Then, 
we  compute  flow  over  a  thin  cylinder  (wire)  and  compare  computed 
and  measured  flow  properties.  Next,  we  compute  flow  over  two 
axisymmetric  bodies  using  both  measured  surface  and  boundary-layer- 
edge  pressure  distributions  and  deduce  that  including  9p/3y  in  a 
thick  boundary  layer  computation  will  account  for  most  of  the  observed 
differences  from  thin  layers.  Finally,  we  devise  an  integral  method 
for  computing  3p/3y  in  a  thick  axisymmetric  boundary  layer. 

5.1  PERTURBATION  ANALYSIS 


The  two  most  noteworthy  scaling  schemes  for  axisymmetric  boundary 
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layers  are  those  of  Rao  J  and  Richmond  .  In  Subsections  5.1.1 
and  5.1.2  we  use  perturbation  methods  to  show  that  the  Wilcox- 
Rubesin  model  is  consistent  with  Rao's  scaling  in  the  viscous  sub¬ 
layer  and  with  Richmond's  scaling  in  the  wall  layer. 


5.1.1  Sublayer  Scaling 


For  an  axisymmetric  boundary  layer  with  constant  pressure,  close 
to  the  surface  the  Wilcox-Rubesin  model  simplifies  to: 


r(v  +  y*e) 


du 

dr 


2 

r  u 

o  T 


(19) 


v  +  °*e)rg.] 

v  +  0E)r£r] 


Y*e  (Jii)2  _  0*u)k  (20) 

dr 

YW  -  [e  +  20  (H7)2]  “3  (21) 


where  v  and  e  denote  molecular  and  eddy  diffusivity,  r  is  radial 
distance,  u  is  velocity,  rQ  is  body  radius  and  uT  is  friction 
velocity . 


Deep  in  the  sublayer,  the  eddy  viscosity  can  be  neglected,  i.e., 
e  <<  v.  Then,  if  we  nondimensionalize  according  to: 


K  *  k/u2 


W  =  vu/u' 


L  =  u^fc/v 


r  *  u  r  /v 

O  TO 

Y  =  r*  log  (r/r  ) 


(22) 


substituting  Equations  (22)  into  Equations  (19-21)  yields 

dU/dY  =1  (23) 

d2K/dY2  =  y*K/W  -  6*WKexp(2Y/r*)  (24) 

d2W2/dY2  =  yW  -  BW3exp(2Y/r*)  -  2a(dL/dY)2W3  (25) 


For  the  most  experimental  measurements  made  on  flows  of  this  type 
the  non-dimensional  radius,  r*,  is  very  large  so  that 

exp(2Y/r*)  =  1  +  0(l/r*)  (26) 

Thus,  in  the  limit  r+  >>  1,  just  as  in  the  two-dimensional  case  it 

is  easily  verified  that  the  production  terms,  y*K/W  and  yW,  and 
the  term  in  Equation(25)  proportional  to  (dL/dY)2  are  negligible. 
The  solution  to  Equations  (23-25)  is  then  trivially  shown  to  be 

U  =  Y  v 

W=  20/( 8Y2 )  [  (27) 

E  ~  Y4  ) 

Writing  the  dimensional  equivalents  of  Equations  (27)  yields  the 
following: 
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v- 


u/uT  =  r^log(r/ro) 

vu>/u2  =  20v/  { 8  ( log  £-)}2 

o 

k/u*  -  <r+log  §-)' 

o 

The  velocity  profile  predicted  in  Equations  (28)  is  precisely  that 
given  by  theRao  scaling.  To  test  the  limiting  u/u^  asymptotic  form, 
we  have  examined  results  of  a  numerical  computation  for  flow  past 
a  thin  cylinder  aligned  axially  with  the  freestream  (complete  detail 
of  the  computation  are  given  in  Subsection  5.2).  For  the  data  sum¬ 
marized  below  in  Table  1  we  have  r+  =  667*8.  As  indicated,  although 

0  + 

the  effect  is  small  because  of  the  large  value  of  r  ,  the  numerical 
u/uT  agrees  with  the  values  predicted  by  the  first  of  Equations  (28) 
to  three  decimal  places.  Thus,  the  analysis  is  self  consistent. 
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Table  1.  Comparison  of  Computed  &  Predicted 
Sublayer  Velocity  Profiles 


uTy/v 

r^log(r/rQ) 

(u/ux  Computed 

0.6056 

0.6053 

0.6054 

1.235 

1.234 

1.234 

1.889 

1.886 

1.886 

2.568 

2.563 

2.563 

3.274 

3.266 

3.266 

4.007 

3-995 

3-995 

4.769 

4.752 

4.752 

5.1.2  Wall-Layer  Scaling 

We  turn  now  to  the  wall  layer,  viz,  the  region  sufficiently  close  to 
the  surface  for  convective  terms  to  remain  negligible  yet  far  enough 
that  the  molecular  viscosity  can  be  neglected  relative  to  the  eddy 
viscosity.  This  is  the  classical  law-of  the  wall  region  for  planar 
boundary  layers.  We  again  nondimensionalize  according  to  Equations 


(22)  with  the  exception  that  our  normal  coordinate  is  redefined  a; 


follows . 


+  +  + 

Y  =  y  (1  +  y  /rQ)  where  y  =  (r  -  r  ) 


(29) 


The  wall  layer  is  most  conveniently  analyzed  by  treating  U  as  the 
dependent  variable.  Dropping  molecular  viscosity  in  Equations 
(19-21)  i  the  equations  for  the  wall  layer  are: 


(1  +  2Y/r^ )  (K/W)dU/dY  =  1 


a*d2K/dU2  =  6*(1  +  2Y/r^)K2  -  1 


(30) 


ad2W2/dU2  =  6(1  +  2Y/r*)Kw2  -  yW2/K  +  2o (dL/dU) 2W Vk  (32) 

In  the  limit  of  very  large  r*  these  equations  are  Identical  to  those 
for  planar  boundary  layers.  From  Equation  (31)  there  follows 
immediately : 


E  =  1//6*  +  0(l/r*) 


(33) 


Substituting  Equation  (33)  into  Equation  (32)*  the  solution  for  W 
becomes  (noting  that  the  closure  coefficients  have  been  chosen  to 
insure  that  y  =  6/6*  -  2ok2//&w  where  <  Is  Karman's  constant): 


W  =  -i 


■c/g* 


exp(-  kU) 


Finally,  substituing  Equations  (33-3^)  into  Equation  (30)  yields 
the  velocity  profile,  viz, 

U  =  7  logY  +  0( 1/r*)  (■ 

or  in  dimensional  form  we  have: 


u/uT  =  ^  logy+  +  . . . 


where 


(37) 


y  =  y  (l  +  y  /2rQ) 

Equations  (36-37)  are  precisely  those  of  the  Richmond  scaling. 

In  Figure  1  we  compare  our  numerical  predictions  for  the  thin  clinder 
computation  of  the  next  subsection  with  Equations  (36-37).  In  the 
comparison  we  have  included  Coles’  wake  component  so  that  the  complete 
velocity  profile  becomes: 

u/uT  =  i  logy+  +  B  +  —  sin2  (75-  ^ )  (38) 

6 

where  we  use  k  =  .Ml,  B  =  5-0  and  tt  =  -.1M.  Note  that  the  inferred 
wake  strength  7  of  -.lM  implies  that  the  transverse  curvature  effect 
is  similar  to  a  favorable  pressure  gradient.  This  point  is  consistent 
with  the  observations  of  Patel-1-?  The  overall  agreement  between  the 
numerical  profile  and  Equation  (38)  verifies  the  perturbation  analysis. 

5.2  FLOW  PAST  A  WIRE 

In  this  subsection  we  present  results  of  a  numerical  computation  for 

flow  past  a  "wire",  viz,  a  thin  cylinder  aligned  axially  with  the 

freestream.  The  case  we  have  chosen  is  that  experimentally  investi- 
16 

gated  by  Yu  in  which  the  Reynolds  number  based  on  cylinder  radius  is 

Rer  =  15,250  (39) 

o 

Unit  Reynolds  number  for  the  flow  is  1.83*105  and  measurements  were 
made  between  axial  distances,  x,  from  2  feet  to  8  feet.  At  x  =  8 
feet  the  ratio  of  boundary  layer  thickness  to  ro  is  approximately  2 
so  that  this  flow  includes  a  moderately  "thick"  boundary  layer. 

Computation  was  initiated  at  x  =  0  from  laminar  starting  conditions. 

The  freestream  turbulence  intensity  was  adjusted  to  yield  a  match 
between  computed  and  measured  momentum  thickness  at  x  =  2  feet. 

Figure  11  compares  predicted  boundary  layer  thickness,  6,  momentum 
thickness,  9,  shape  factor,  H,  and  skin  friction,  cf,  with  corres¬ 
ponding  measurements.  The  data  shown  for  6,  6,  and  H  are  the  original 
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Figure  10.  Comparison  of  a  numerical  velocity  profile  with  that  predicted  by 
perturbation  analysis  of  the  Wilcox-Rubesin  model. 


COMPUTED 

PLANAR  BOUNDARY  LAYER 


Figure  11.  Comparison  of  computed  and  measured  integral  properties  for  flow  past  a  thin 
cylinder  aliened  axially  with  the  frcestream;  Re^  =  1  '3 , 2 f  0 . 


data  reported  by  Yu.  The  skin  friction  data  have  been  obtained  from 
Clauser  plots  rather  than  using  Yu*s  Preston  tube  measurements  which 
are  thought  to  be  in  error1^. 

As  shown,  computed  and  measured  shape  factor  and  skin  friction  differ 
by  less  than  5%.  Computed  boundary-layer  thickness  is  about  10% 
lower  than  measured  while  the  numerical  momentum  thickness  ultimately 
is  about  15%  lower  than  measured.  The  fact  that  our  skin  friction 
is  so  close  to  the  measured  values  while  the  momentum  thickness 
discrepancies  are  much  larger  leaves  us  with  the  same  concern  expres¬ 
sed  by  Patel1^ regarding  the  two-dimensionality  of  the  experimental 
flowfield.  The  figure  also  shows  computed  results  for  a  corresponding 
planar  boundary  layer.  As  indicated,  the  primary  effects  of  transverse 
curvature  are  to  increase  skin  friction  and  momentum  thickness  and  to 
decrease  boundary-layer  thickness  and  shape  factor. 

Figure  12  compares  the  computed  and  measured  velocity  profiles  at  the 
axial  location  x  =  8  feet.  With  exception  of  the  points  lying  between 
y+  =  150  and  400,  the  computed  velocity  profile  fares  smoothly  through 
the  data.  Even  for  the  four  data  points  lying  between  y+  =  150  and 
400,  the  maximum  difference  is  less  than  5%  •  Inspection  of  the 
profile  for  the  corresponding  planar  case  shows  that  the  effect  of 
transverse  curvature  is  to  alter  the  slope  of  the  profile,  a  trend 
consistent  with  that  displayed  in  the  data. 

Thus,  for  this  constant-pressure  flow,  the  model  equations  pre¬ 
dict  flow  properties  which  are  reasonably  close  to  corresponding 
measurements.  As  will  be  seen  in  the  next  subsection,  thick 
boundary  layers  experiencing  an  adverse  pressure  gradient  cannot 
be  as  accurately  computed  if  we  restrict  ourselves  to  classical 
thin-shear-layer  approximations. 


Figure  12.  Comparison  of  computed  and  measured  velocity  profiles  for  flow  past  a  thin 
cylinder  aliened  axially  with  the  freestream;  Re  =  lF,2r>0,  x  =  8  feet. 


5.3  BODIES  WITH  ADVERSE  PRESSURE  GRADIENT 


In  this  Subsection,  we  make  a  "first  cut"  at  two  thick  axisymmetric 
boundary  layers  with  pressure  gradient  and  streamline  curvature. 

We  emphasize  that  our  computations  are  indeed  only  a  first  cut  as 
experimental  measurements  indicate  nontrivial  variation  in  pressure 
across  the  boundary  layers  while  we  have  used  a  standard  boundary- 
layer  program  which  assumes  9p/8y  =0.  In  order  to  provide  a  first 
estimate  of  the  importance  of  having  variable  pressure  across  the 
layer,  we  have  done  both  computations  first  using  the  measured  edge- 
pressure  distribution  and  then  with  the  measured  surface-pressure 
distribution . 

18 

The  two  cases  considered  are  the  "modified  spheroid"  and  the 
19 

"low-drag  body"  experimentally  evaluated  by  Patel.  In  Subsection 
5.3.1  we  compare  computed  and  measured  integral  properties.  Sub¬ 
section  5-3.2  compares  computed  and  measured  velocity  and  Reynolds 
shear-stress  profiles  for  the  low-drag  body.  Next,  in  Subsection 
5.3.3,  we  compare  predicted  peak  eddy  viscosity  and  mixing  length 
with  corresponding  measurements.  Finally,  Subsection  5-3.^  examines 
the  predicted  effect  of  streamline  curvature  for  the  modified 
spheroid . 

5.3.1  Integral  Properties 

Figures  13  and  14  compare  computed  and  measured  integral  properties 
for  the  two  axisymmetric  bodies.  Note  that  the  modified  spheriod 
computation  has  been  initiated  at  x/L  =  .662  using  the  experimentally 
measured  profiles  to  establish  starting  profiles.  Reynolds  number 
based  on  body  length  is  Re^  =  1.26*10^.  For  the  modified  spheriod, 
boundary-layer  thickness,  6,  exceeds  the  body  radius,  r  ,  for  all 
points  beyond  about  x/L  =  .925.  Similarly  for  the  low-drag  body, 
computation  has  been  Initiated  at  x/L  =  .6,  again  using  the  experimen¬ 
tally  measured  profiles  for  starting  conditions.  Reynolds  number, 

Re^,  is  1*20*10^  while  6  exceeds  rQ  beyond  about  x/L  =  .85. 

Focusing  first  on  the  modified  spheroid,  note  that  the  boundary  layer 
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Fipure  .  Comparison  of  computed  and  measured  integral  properties  for  Patel 
low-dran  body;  ReT  =  1.20'106. 


remains  attached  to  x/L  =  .99  when  we  use  Cpg  and  separates  at  x/L  = 

. 992  when  we  use  Cpw«  For  both  pressure  distributions,  the  predicted 
growth  of  the  boundary-layer  thickness  is  reasonably  close  to  the 
measured  trend  with  the  experimental  data  beyond  x/L  =  .9  lying  between 
the  two  computed  curves.  Similarly,  the  experimental  skin-friction, 

Cf,  data  tend  tc  fall  between  the  CPe  and  Cpw  curves,  particularly  in 
the  region  beyond  x/L  =  .9  where  9p/9y  is  known  from  measurements  to 
be  nontrivial.  On  balance,  the  Cpw  computation  predicts  momentum- 
thickness,  0,  and  shape-factor,  H,  distributions  which  are  closer  to 
the  corresponding  measurements  in  the  interesting  (i.e.,  x/L>.9) 
portion  of  the  flow. 

Turning  now  to  Patel's  low-drag  body,  note  that  the  experiment  was 
designed  to  avoid  any  complications  attending  incipient  separation, 
such  as  that  observed  for  the  modified  spheroid.  As  shown  in  Figure 
19  neither  computation  predicts  separation  and  differences  between  the 
two  computations  are  less  dramatic  than  in  the  modified  spheroid  case. 

As  with  the  preceding  computation,  predicted  overall  growth  of  the 
boundary-layer  thickness  is  quite  close  to  the  measured  growth.  For 
both  computations,  computed  and  measured  Cp  differ  by  less  than  1251. 

As  with  the  modified  spheroid,  differences  between  computed  and  measured 
6  and  H  distributions  are  smallest  for  the  Cp  computation. 

In  summary,  comparison  of  predicted  and  measured  integral  properties 
indicates  that  the  experimental  6  and  Cp  fall  between  the  limiting  Cpe 
and  CPw  computations.  Predicted  0  and  H  distributions  are  closest  to 
the  corresponding  measurements  when  we  use  Cp  =  Cp  .  Hence,  it  appears 
that  in  taking  a  closer  look  at  our  predictions,  the  Cr.  computations 
will  provide  more  insight  than  those  using  Cpe- 

5.3.2  Velocity  and  Reynolds-Stress  Profiles 

Because  the  modified  spheroid  computation  separates,  there  is  little 
positive  information  to  be  gleaned  from  detailed  comparison  of  computed 
and  measured  profiles.  Hence,  in  this  Subsection,  we  focus  our  atten¬ 
tion  on  the  low-drag  body. 


Figure  15  compares  computed  and  measured  velocity  profiles  at  four  axial 
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Figure  15.  Comparison  of  computed  and  measured  velocity  profiles  for  Patel’s 
low-drag  body;  ReT  =  1.20*106. 


stations,  viz,  x/L  =  .80,  .88,  .92,  .96.  As  shown,  the  computed 
boundary  layer  has  experienced  stronger  deceleration  than  in  the 
experimental  case.  This  is  unsurprising  as  we  have  used  the  surface- 
pressure  distribution  throughout  the  layer  which  is  attended  by  a 
larger  gradient.  Clearly,  the  overall  growth  matches  that  measured. 

Even  at  x/L  =  .9 6,  differences  between  computed  and  measured  velocities 
are  less  than  7%  of  scale. 

Figure  16  compares  computed  and  measured  Reynolds  shear-stress  profiles. 
Discrepancies  between  the  profiles  at  x/L  =  .80  are  surprising,  parti¬ 
cularly  in  the  light  of  the  close  agreement  between  the  velocity 
profiles  at  this  station.  Overall,  the  numerical  shear-stress  profiles 
are  within  about  15%  of  scale  of  the  measured  profiles.  The  shapes 
and  fullness  of  the  experimental  profiles  are  reasonably  well  simulated 
throughout  the  tail  region  of  the  body. 

5.3.3  Turbulence  Properties 

In  this  subsection  we  concentrate  on  two  key  features  of  the  turbulence 
field  which  have  not  been  accurately  simulated  with  simpler  turbulence 
models.  Specifically,  we  examine  the  variation  of  peak  mixing-length, 
lmax’  and  edcJy  viscosity,  Figure  17  compares  our  Cpw  predictions 

with  values  inferred  from  Patel's  data.  As  shown,  for  both  bodies, 
model-predicted  l  and  e  fall  off  rapidly  as  we  approach  the  tail 
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of  the  body.  Considering  the  inaccuracies  attending  differentiation 

of  the  experimental  data  required  to  infer  A  and  e__  .  our  predic- 

^  max  max  ^ 

tions  must  be  considered  well  within  the  error  band  of  the  data. 

5.3-^  Streamline-Curvature  Effects 

To  obtain  an  estimate  of  the  effects  of  streamline  curvature  on  our 

predictions,  we  have  rerun  the  Cpw  modif ied-spheroid  computation 

without  account  for  streamline  curvature.  All  results  presented  in 
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Subsections  5* 3.1-5. 3* 3  Include  the  Wilcox-Chambers  curvature  modi¬ 
fication  to  the  turbulent  mixing-energy  equation.  Figure  18  shows  effect 
on  the  velocity  profile  at  x/L  =  .93.  For  the  computation  with  no 
curvature,  the  numerical  profile  passes  through  all  of  the  measured 
points.  With  curvature  included,  the  velocity  profile  shows  additional 
deceleration,  consistent  with  the  effect  of  convex  curvature.  Note 


WITH  CURVATURE 
NO  CURVATURE 


PATEL 


that  the  close  agreement  between  results  of  the  no-curvature  computa¬ 
tion  and  the  experimental  data  is  fortuitous  as  the  numerical  pressure 
gradient  is  larger  through  much  of  the  boundary  layer  than  it  is  in 
the  experimental  case.  Overall,  the  effect  of  curvature  for  the 
modified  spheroid  is  small,  giving  rise  to  changes  in  integal  pro¬ 
perties  of  less  than  5% 

In  summary,  results  of  these  two  computations  indicate  the 

following. 

1.  With  no  modifications  to  the  Wilcox-Rubesin 
turbulence  model,  computed  boundary-layer 
thickness  virtually  duplicates  measured 
thickness  for  both  bodies. 

2.  Performing  the  computations  with  first  the 
measured  boundary-layer-edge  pressure 
distribution  and  then  with  the  measured 
surface  pressure  distribution  (neglecting 
pressure  variation  across  the  layer  in  each 
computation)  yields  skin-friction  distribu¬ 
tions  which  overall  fall  above  and  below 
measured  skin  friction,  respectively. 

3.  Shape  factor  and  momentum  thickness  are  in 
closer  agreement  with  corresponding  measure¬ 
ments  when  the  wall  pressure  is  used. 

4 .  For  one  of  the  two  bodies,  our  predictions 
suggest  that  streamline  curvature  plays  a 
relatively  minor  role. 

Thus,  we  conclude  that  the  primary  differences  between  the  "thick" 
and  the  "thin"  boundary  layer  are  primarily  caused  by  the  pressure 
gradient  across  the  former,  i.e.,  3p/3y.  In  the  next  Subsection, 
we  devise  a  straightforward  method  for  including  3p/3y  in  a  standard 
boundary-layer  computation. 


5.M  PRESSURE  VARIATION  IN  A  THICK  BOUNDARY  LAYER 

Results  of  the  preceding  Subsection  demonstrate  the  importance  of 
pressure  variations  across  a  thick  axisymmetric  boundary  layer. 
Recapitulating,  we  found  that  in  using  a  conventional  boundary 


layer  program,  results  obtained  differ  depending  upon  the  pressure 
distribution  used.  On  the  one  hand,  using  the  measured  freestream 
pressure  distribution  tends  to  underpredict  momentum  thickness  and 
to  overpredict  skin  friction.  On  the  other  hand,  using  the  measur¬ 
ed  surface  pressure  tends  to  reverse  the  situation  with  momentum 
thickness  overpredicted  and  skin  friction  underpredicted.  Pre¬ 
sumably,  permitting  the  pressure  to  vary  across  the  layer  will 
yield  numerical  predictions  which  fall  somewhere  between  these 
two  limits.  Because  the  experimental  data  lie  between  these  two 
limits,  predictions  should  then  lie  quite  close  to  corresponding 
measurements.  The  object  of  this  Subsection  is  to  devise  a  method, 
compatible  with  a  conventional  boundary-layer  marching  algorithm, 
for  predicting  the  pressure  variation  through  a  thick  boundary  layer 

In  order  to  account  for  nonzero  3p/3y  in  a  boundary  layer,  we  must, 

in  principle,  solve  the  vertical  as  well  as  the  streamwise  momentum 

equation.  For  a  thick  boundary  layer  on  a  surface  with  curvature 
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k,  Patel  has  shown  that  the  vertical  momentum  and  continuity 
equations  are: 

u  3v/3x  +  hv3v/3y  -  ku2  +  ^  3p/3y  =  0  (40 

3(ru)/3x  +  3(rhv)/3y  =0  (41 

where  u  and  v  denote  streamwise  and  vertical  velocity  components, 
x  and  y  are  streamwise  and  normal  distances, p  is  density,  p  is 
pressure,  r  is  radial  distance,  and  h  is  the  metric  defined  by 

h  =  1  +  icy  (42 

The  curvature,  k,  is  understood  to  be  positive  for  convex  surfaces 
and  negative  for  concave  surfaces. 

One  approach  to  solving  for  the  pressure  would  be  to  discretize 
Equation  (40)  and  solve  for  p  once  v  is  known  from  the  standard 
boundary-layer  solution.  This  is  not  very  convenient  when  the  pro¬ 
gram  is  formulated  in  terms  of  Levy-Lees  variables,  however,  and 
all  of  our  programs  use  these  variables.  As  an  alternative,  we  have 


chosen  to  implement  an  integral  solution  approach.  Proceeding 
in  the  classical  manner,  we  define  the  following  boundary-layer 
thickness  parameters: 


Displacement  Thickness 
U-Momentum  Thickness.. 
V-Momentum  Thickness.. 


r  6* 

-'S 

(1-u/u  )rdy 

(43) 

r  0 

-'S 

~  (1-u/u  )rdy 

0 

(44) 

r6 

V 

-  'o 

~  (l-v/vg)rdy 

(45) 

where  6  is  the  boundary-layer  thickness.  Then,  introducing  the 
following  quantities: 


r  =  r  +  6cos$  ,  h  =  1  +  <6,  r  6  =  f  rdy 

60  9  6  *or  o  ° 


the  integral  form  of  Equations  (40-41)  is  as  follows. 


re  theve  ~  ued(S/dx3  +  d[ro(6r  "  6#)  ue  ^  /dx  =  0 


d  [r  u  v  0  ]/dx  -  r  (5  -  <$*)u  dv  /dx  +  ter  u2(6  -6*-0) 

l  o  e  6  o  x*  6  6  o  6  r 


=  f6o  £(9p/3y)rdy 


In  order  to  make  further  headway  we  must  postulate  profiles  for  v 


and  p.  Close  examination  of  data  for  two  thick  axisymmetric  bodies 
shows  that  at  all  stations  for  which  data  were  tabulated,  both  v 
and  p  can  be  closely  approximated  with  simple  linear  profiles,  viz, 
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v  =  ve(y/o)  and  p  =  pg  +  (pw~Pe )  (l-y/6 ) 


Thus,  the  definition  of  0^  is  altered  in  the  obvious  way  while 


Equation  (48)  simplifies  to: 


d [r  u  v  0  ]/dx  -  r  (6  -6*)u  dv  /dx  +  <r  u2(6  -6#-9) 

1-0  6  6  O  i*  66  061? 

_  /  t:2  r  \  T  r  i  2  .  f  /  p  f  /  h  \  1  (  r> —  r\ \ 


(rou:/26)[6rt|<6(6r-6/4)]  (Cp£  -Cpw) 


In  order  to  test  this  formulation,  we  have  applied  Equations  (47 
and  50)  to  three  flows,  in  each  case  working  with  measured  flow 
properties.  The  three  cases  are  Yu's1®  flow  over  a  wire,  Patel's 
modified  spheroid1®  and  Patel's  low-drag  body1^,  respectively. 

For  Yu’s  wire  flow,  we  find  that  the  difference  between  edge  and 
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surface  pressure  coefficients  is  of  order  10  ,  as  would  be  ex¬ 

pected  for  this  constant  pressure  case.  Figures  19  and  20  compare 
predicted  and  measured  pressure  coefficient  difference  for  the 
two  Patel  bodies.  Considering  the  numerical  crudeness  attending 
differentiation  of  the  data,  particularly  v  ,  the  agreement 
between  the  predictions  of  Equation  (50)  and  the  measured  C 

P 

difference  is  excellent. 


6.  DISCUSSION 
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We  have  developed  a  new  three-dimensional  boundary-layer  program, 
EDDY3,  suitable  for  application  to  arbitrary  ship  hulls.  The  pro¬ 
gram  embodies  the  Wilcox-Rubesin  two-equation  turbulence  model  and 
uses  an  accurate,  efficient  numerical  procedure  based  on  the  Elot- 
tner  "variable-grid"  method  coupled  with  the  Krause  explicit  mar¬ 
ching  algorithm. 

Results  described  in  Sections  3  and  A  show  that  very  accurate  num¬ 
erical  results  can  be  obtained  with  relatively  coarse  finite-dif¬ 
ference  grids.  Computing  times  are  quite  modest  as  a  complete  ship 
hull  computation  can  be  done  with  about  a  half-hour  of  UN'IVAC  1108 
execution  time.  This  corresponds  to  approximately  three  minutes  on 
a  CDC  7600  computer.  Additionally,  the  numerical  algorithm  dis¬ 
plays  no  noticible  difficulty  when  reverse  crossflow  is  present. 

Comparison  of  computed  and  measured  flow  properties  for  the  SSFA 
Model  720  and  the  HSVA  Tanker  show  that,  on  the  one  hand,  our  pre¬ 
dictions  are  quite  a  bit  closer  to  measurements  than  these  reported 
in  the  1980  SSPA-ITTC  Workshop  on  Ship  Boundary  Layers.  On  the 
other  hand,  differences  are  too  large,  particularly  near  the  ship 
stern,  for  us  to  claim  the  problem  is  solved.  The  fact  that  the 
boundary  layer  becomes  "thick"  approaching  the  stern  is  no  doubt 
the  cause  of  the  differences  observed  between  theory  and  experiment 

Our  analysis  of  thick  axisymmetric  boundary  layers  in  Section  5 
shows  that,  in  order  to  obtain  accurate  "thick"  boundary-layer  pre¬ 
dictions,  accounting  for  the  normal  pressure  gradient,  3  p / 3  y  , 
probably  is  all  that  is  needed  above  and  beyond  conventional  thin- 
shear-layer  approximations.  Streamline  curvature  variation  appears 
to  play  a  relatively  minor  role  in  the  "thick"  boundary  layer. 

Additional  research  is  needed  to  confirm  the  importance  of  3p/3y  in 
ship-hull  computations.  The  formalism  developed  in  Subsection  5.“ 
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can  be  easily  generalized  for  a  3-D  boundary  layer  and  should  be 
included  in  EDDY3 .  We  could  then  repeat  the  two  ship-hull  com¬ 
putations  of  Section  to  confirm  our  hypothesis.  Additionally, 
it  would  be  instructive  to  use  measured  freestream  flow  condition 
in  the  computations. 
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It  contains  complete  program  listings  and  input/output  descrip¬ 
tions  for  the  three  programs  SHPMSH,  VELOC  and  EDDY3.  To  obtain 
a  copy  of  this  report  including  the  Appendix,  contact  either  the 
David  W.  Taylor  Naval  Ship  Research  and  Development  Center  or 
DCW  Industries,  Inc. 
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